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Abstract 

We study the dynamical and statistical behavior of the Hamiltonian Mean Field 
(HMF) model in order to investigate the relation between microscopic chaos and 
phase transitions. HMF is a simple toy model of fully-coupled rotators which 
shows a second order phase transition. The canonical thermodynamical solution is 
briefly recalled and its predictions are tested numerically at finite A^. The Vlasov 
stationary solution is shown to give the same consistency equation of the canonical 
solution and its predictions for rotator angle and momenta distribution functions 
agree very well with numerical simulations. A link is established between the behav- 
ior of the maximal Lyapunov exponent and that of thermodynamical fluctuations, 
expressed by kinetic energy fluctuations or specific heat. The extensivity of chaos in 
the — > oo limit is tested through the scaling properties of Lyapunov spectra and 
of the Kolmogorov-Sinai entropy. Chaotic dynamics provides the mixing property 
in phase space necessary for obtaining equilibration; however, the relaxation time to 
equilibrium grows with A'^, at least near the critical point. Our results constitute an 
interesting bridge between Hamiltonian chaos in many degrees of freedom systems 
and equilibrium thermodynamics. 
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1 Introduction 



Many-particle systems can show collective behavior when the average kinetic 
energy is small enough. This collective macroscopic behavior can coexist with 
chaos at the microscopic level. Such a behavior is particularly evident for sys- 
tems that have a phase transition, for which a nonvanishing order parameter 
measures the degree of macroscopic organization, while at the microscopic 
level chaotic motion is a source of disorder. The latter can induce non triv- 
ial time dependence in the macroscopic quantities, and it would be desirable 
to relate the time behavior of such quantities and their fluctuations to the 
chaotic properties of microscopic motion, measured through the Lyapunov 
spectrum. A naive idea is that an increase of chaos as the energy (tempera- 
ture) is increased should be accompanied with a growth of fluctuations of some 
macroscopic quantity. These should be maximal at the critical point and then 
drop again at high energy. In this paper we study a model of fully-coupled 
Hamiltonian rotators which realizes such a behavior, it has been called Hamil- 
tonian Mean Field (HMF) model [1,2]. It can also be considered as a system 
of interacting particles moving on a circle. This system has a second order 
phase transition and in the ordered phase the rotators are clustered; the high 
temperature phase is a gaseous one, with the particles uniformly distributed 
on the circle. It has been shown in ref. [2] that the maximal Lyapunov ex- 
ponent grows up to the critical energy density Uc and then drop to zero in 
the whole high temperature phase in the N ^ oo limit. Correspondingly one 
observes a growth of kinetic energy fluctuations up to the critical point and 
then a phase of vanishing fluctuations. Finite effects complicate this sim- 
ple picture. In the high temperature phase the maximal Lyapunov exponent 
vanishes quite slowly (with A"~^/^) and finite size effects influence the first re- 
gion below the critical point. In this region the system displays metastability: 
starting far from equilibrium, this is reached in a time Tr which grows with A^. 
On the contrary, the extremely low energy phase is characterized by a weak 
N dependence, with the maximal Lyapunov exponent Ai which behaves as 
Ai ~ VU. 

Although the model is extremely simplified, it shares many features with more 
complex models, for which the relation between chaotic motion at the micro- 
scopic level and collective macroscopic properties has been studied. Let us 
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mention studies in solid state physics and lattice field theory [3-7]. However, 
it has been actually in nuclear physics [8,9], where there is presently a lively 
debate on multifragmentation phase transition [8-16], that the interest in the 
connection between chaos and phase transitions has been revived. In this case 
in fact, an energy /temperature relation quite close to the HMF model has been 
observed [10] and critical exponents have been measured experimentally [11]. 
Statistical thermodynamical models [12] and percolation approaches [14] have 
been proved to give a good description of the experimental data, though the 
dynamics is missing. On the other hand classical molecular dynamics models 
[9,13,15] seem to contain all the main ingredients, but have the disadvantage 
that a detailed understanding of the dynamics can be too complicated. In this 
respect, the HMF model can be very useful in clarifying some general dynam- 
ical features which could be eventually compared with real experimental data. 
In fact, when studying nuclear multifragmentation, one deals with excited 
clusters of 100-200 particles interacting via long-range (nuclear and Coulomb) 
forces. Quantum effects are relevant only at very low energy. In fact in the 
nuclear case, at very low energy, T is not linear in U , but grows as \/U be- 
cause nucleons are fermions [10,12]. However, a classical picture should be 
quite realistic in the critical region where the excitation energy is substantial 
[16]. 

In this paper we present new numerical data concerning both statistical quan- 
tities, like specific heat and distribution functions, and chaotic probes, like 
Lyapunov spectra and Kolmogorov-Sinai entropy. Moreover, we add to the 
theoretical analysis of the model a thorough treatment of differences in the 
fluctuating quantities between the canonical and microcanonical ensembles. 
We also investigate in detail the relaxation to equilibrium and compare numer- 
ical results with a complete self-consistent Vlasov calculation of distribution 
functions. Finally a comparison of numerically obtained maximal Lyapunov 
exponents with theoretical formulas is attempted. 

The paper is organized as follows. In Sec. 2 we briefly discuss the details of the 
HMF model. The equilibrium statistical mechanics and the continuum Vlasov 
solution are described in Sects. 3 and 4 respectively. In Sec. 5 we discuss the 
relaxation to equilibrium and in Sec. 6 we present the numerical calculations 
of the Lyapunov spectra and Kolmogorov-Sinai entropy as a function of the 
energy and A^. Analytical estimates are discussed in Sec. 7 and conclusions 
are drawn in Sec. 8. 
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2 The HMF model 



The Hamiltonian we consider is the following 

H{{ei},{p,})^K + V , (1) 



where 



i — 1 ijj — 1 



are the kinetic and potential energies. The model describes the motion of N 
rotators characterized by the angle 9i G [0, 27r): each rotator interacts with all 
the others. One can define a spin vector associated with each rotator iiij ~ 
{cos{9i), sin{9i)) . The Hamiltonian then describes N classical spins similarly to 
the XY model, and a ferromagnetic or an antiferromagnetic behavior according 
to the positive or negative sign of e [1]. In the following we will consider only 
the ferromagnetic (attractive) case and moreover we put e = 1 without loss of 
generality. After defining M = Th=\ ^ = {Mx, My) , the potential energy 
V can be rewritten as 

^ = y (1 - (M'x + Ml)) = ^(1 - M^) . (3) 



The equations of motion then read 

Oi^Pi , Pi^ -sin(ei)Mx + cos(ei)My , i ^ 1, N . (4) 

The evolution equations of the tangent vector are 

Se, = Spi , 6p, = -^-^^6e, , t = l,...,N. (5) 

j i j 



where the diagonal and off-diagonal terms of the Jacobian = —d'^V/dOidO 



are 
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Jii = -cos{ei)Mx - sin{ei)My + (6) 
Jij = ^cos{0i - 0j) , iy^j . (7) 
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Expression (6) can also be written for convenience as: 

Jii = -McosiOi - 0) + 4 



(8) 



where is the phase of M . 

The HMF model was initially proposed in Ref [1]. A time discrete version 
of it was previously introduced in Ref. [18] in the form of globally coupled 
Chirikov standard maps. The Lyapunov instabihty of the HMF has been first 
studied numerically in Ref. [6], and analytically in Ref. [19] using a Riema- 
niann geometry approach [7]. The connection between Lyapunov instability 
and thermodynamical properties has been established in Ref. [2]. In the orig- 
inal proposal of the HMF [1], the relation of the model with self-gravitating 
systems (in the attractive case) and charged sheets models (in the repulsive 
case) was studied; this has been also taken over for the gravitational case in 
Ref. [20] and in the context of plasma models [21,22]. It has also been shown 
that the HMF is the peculiar representative of a larger class of models, whose 
thermodynamics can be solved exactly [21]. Moreover, anomalous diffusion 
properties of a generalization of the model, with two angles for each rotator, 
have been recently studied [23] . 



3 Canonical and microcanonical results 

It is interesting to look at the predictions of statistical mechanics. We restrict 
here to the N ^ oo limit, where microcanonical and canonical ensemble results 
coincide for averaged quantities (apart from exceptions near first order phase 
transitions, see Rcfs. [24]). The free energy of our model can be easily obtained 
in the canonical ensemble. This was done in Ref. [1] and the result reads 



where P — l/ksT (the Boltzmann constant ks is set to 1) and is the 

modified Bessel function of i-th order. The auxiliary variable y is introduced 
to decouple the particles by the usual Hubbard-Stratonovich trick for Gaussian 
integrals, and the search of the maximum in Eq. (9) (which is a consequence 
of the continuum limit — > oo solved with the saddle point method) leads 
to the consistency equation 



1 /27r\ 3 f 



+ log(27r/o(|/))) , (9) 




(10) 
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The magnetization M — Ii/Iq is obtained by solving the consistency equation, 
which is done numerically. It vanishes for j3 < j3c = 2, the inverse critical 
temperature. At low temperatures it approaches 1 (the limit of I\/ Iq). 



The energy-temperature relation (sometimes called the caloric curve ) is 



Close to the critical point /3 — > magnetization and energy behave as 



Hence, M vanishes with the 1/2 classical mean field exponent and the specific 
heat Cy = dU/dT is finite at the critical point Cv{Pc) = 5/2 and constant 
(Cy = 1/2) in the high temperature phase. Thus for what concerns the critical 
behavior, Cy ~ (Tc - T)"" with a = 0. 

These theoretical results are compared with those of numerical simulations 
in Fig. 1, where care is taken to use almost equilibrated initial data in order 

to reduce the relaxation time (as discussed in Sects. 4 and 5). We have in- 
tegrated Eqs. (4) using fourth order symplectic algorithms [25] with a time 
step At ~ 0.2, adjusted to maintain the error in energy conservation below 
^ = 10~^. The agreement is good over the whole energy range and finite N 
effects, although present, are weak. 

In the inset of Fig. lb we plot the results obtained starting from non equi- 
librated initial data (the "water bag" initial condition discussed in Sect. 5) 
and although the integration time was quite long (O(IO^)), a sharp disagree- 
ment is observed just below the critical point. A region of negative specific 
heat is present and a continuation of the high temperature phase (linear T 
vs. U relation) into the low temperature one (mctastability). It is very in- 
triguing that these out-of-equilibrium quasi-stationary states (QSS) show a 
caloric curve very similar to the one found for first order phase transitions 
[23,24,26,27]. In that case, however, the corresponding states are equihbrated 
and do not die asymptotically, as it is shown in Sec. 5 for our QSS. In our 
case, at equilibrium, we do not have phase-coexistence, all rotators belong to 
a single cluster (see Sec. 4); phase-coexistence can arise only as a finite N 
non-equilibrium effect. The existence of long- living non-equilibrium states has 




(11) 



thus Uc = Ec/N = 3/4. 




(12) 



(13) 
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been noticed already in Refs. [2,6] and a connection to critical slowing down 
has also been proposed. More recently it has also been found numerically in 
other one- dimensional models [28] and in self-gravitating systems [29], but in 
this case it is not associated to the closeness of a phase transition. The coex- 
istence of different states in the continuum limit near the critical point is a 
purely microcanonical effect. It arises because of the inversion of the t ^ oo 
limit with the N ^ oo one. 



Concerning fluctuations, it is well known that the predictions of microcanoni- 
cal ensemble differ from those of the canonical [30]. For instance, while in the 
canonical ensemble kinetic energy fluctuations are given by 



aUK) _ l {<K^>can-<K>U 1 .... 



their expression in the microcanonical ensemble is [19] 



E = ^ = ^ /i ^ fi5) 

^ Vn l-2M{dM/dT,) 



where = 2 < K > ^ /N (for a rigorous definition of microcanonical temper- 
ature see [31]). The latter is compared with numerical simulations in Fig. 5b. 
The specific heat is [30] 



2 V 



(16) 



and also compares quite well with numerical simulations (see Fig. 5c); although 
finite N effects are obviously stronger for fluctuations than for the averages. 



4 The Vlasov solution 



The Vlasov equation for our system reads 



df df dVdf ^ 
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where f{9,p,t) is the normahzed distribution function and the potential V 
satisfies the equation 



Q2y 2. ~ 



J cos{e-e')f{e',p',t)de'dp' . (is) 



09^ 

-oo 



A simple hypothesis that can be made on the distribution function / is that 
it factorizes as 

f^Mp)9(0,t). (19) 



Although simple, this is a very strong hypothesis and its vahdity can at present 

be justified only by the correctness of the results which are derived thereby 
(see the following). In the high temperature phase one expects a uniform dis- 
tribution in 9 and then / should depend only on p, this fact is then consistent 
with the factorization hypothesis. Lowering the temperature below the tran- 
sition point, a modulation in 9 appears. This is a manifestation of the Jeans 
instability (see Ref. [20]). If the modulation is not too strong, the factorization 
hypothesis is again reasonable. However, in the low temperature/energy phase 
this hypothesis has no clear a priori justification. 



By further requiring that /o is Gaussian 



one gets the equation 



^ + P%-f(^y cos ^ - ^) 9(0) = (21) 



where = J cos 9gd9, My = J sin 9gd9. 

Restricting to the stationary solution, dg/dt — one can easily solve for g{9) 

(22) 



g = go exp 



{My sin 9 + M^ cos 9) 



go exp 



M 

— cos{9 - 0) 



where go = l/(27r/o(M/T)) as imposed by normalization and is the phase 
of M. It must be observed that g is expressed in terms of (M^., My) or (M, (p), 
which are themselves functions of g; so one must solve the problem self- 
consistently, after writing the equations for the two components of the mag- 
netization 
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These equations coincide exactly with the consistency equations (10) of the 
solution in the canonical ensemble. Once (M^., My) are determined by solv- 
ing (numerically) Eqs. (23,24), they can be substituted back into Eq.(22), 
thus obtaining the stationary distribution function. Due to the global phase 
translation invariance of the model, for any M the choice of the phase is ar- 
bitrary, as reflected in the solutions of Eqs. (23,24). With respect to the results 
in Ref. [20] we have fully determined the distribution function by imposing 
the self- consistency conditions through Eqs. (24). 

In the high temperature phase M = thus g — 1/ (27r) is uniform. 

At very low temperature one can use the asymptotic development of Iq 

(25) 



V2 



nz 



to get the Gaussian 



where a"^ = T/M is the variance, which vanishes with T, giving a Dirac-5 at 
zero temperature. 

A comparison of formulas (20,22) with numerical data is shown in Fig. 2; the 
theoretical curve fits the data very accurately with no free parameter. Both 
in the low energy region (Fig. 2a,b) and at higher energy, where the cluster 
drifts (Fig. 2c,d), the agreement is very good. The theory does not determine 
the value of 0, which remains arbitrary; so we have adjusted this value to the 
center of the cluster, which is moving in time quite irregularly. 



5 Slow relcixation to equilibrium 



Around the critical energy, relaxation to equilibrium depends in a very sen- 
sitive way on the adopted initial condition [32]. When starting with "water 
bag" initial conditions, i.e. a fiat distribution function of finite width centered 
around zero for fo{p), and putting all rotators at = {g{9) = 5(0)), we reveal 
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the presence of quasi-stationary (long living) non-equilibrium states (QSS) (see 
the inset of Fig. lb). In Fig. 3 the evolution to equilibrium of the QSS state 
is shown by the time evolution of the reduced distribution function fo{p,t). 
It is only at t — 5 • 10^ that a good reproduction of the Gaussian distribution 
and a convergence to the predicted equilibrium temperature 7)j = 0.4757 for 
U — 0.69 is obtained for = 1000. Such a slow relaxation is observed in the 
region just below the critical point (see Fig. lb) and also around U ~ 1. In 
order to study the A^-dependence of the relaxation time we have roughly 
quantified the distance from the equilibrium state by measuring 

AS = \S{t) -S"'] , (27) 

where S = — J fo{p, t) In fo{p, t)dp is the Boltzmann entropy of the momentum 
reduced distribution function, and S*^^ its equilibrium value when the distribu- 
tion is a Gaussian at the given equilibrium temperature T^. The results, shown 
in Fig. 4, clearly indicate an increase of the relaxation time with N (Fig. 4(a)), 
which can be approximately fitted with a linear law Tr ~ A^. The convergence 
to the equilibrium value of S is not exact at finite and the error decreases 
as N~'^^'^; this is shown in Fig. 4(a) by the convergence to a decreasing value 
(horizontal lines) as A^ is increased from 1000 to 10000. A similar law was 
found in Ref. [1] by studying the time needed to reach equipartition of rota- 
tors velocity in the large U region and in Rcf. [22] analysing the time needed 
to absorb holes in the momentum distribution (so called Dupree structures in 
plasma physics) for the antiferromagnetic HMF. It then seems that various 
indicators agree in suggesting a diverging time scale with A^. However, the 
time scale also depends on U, and what we have here found is that it is much 
greater in the region near the critical point. This could be a manifestation of 
critical slowing down. 



6 Lyapunov spectra and Kolmogorov-Sinai Entropy 

We have computed the Lyapunov spectrum Xi ,i = 1, . . . , N hj the standard 
method of Ref. [33] . The average number of time steps in order to get a good 
convergence was of the order 10^. We discuss in the following numerical results 
for system sizes in between N=10 and N=20000. 

In Fig. 5(a) we plot Ai as a function of U for various A^ values. As the system 
is integrable both in the limit of very small and very large energies (reducing 
in the former case to weakly coupled harmonic oscillators and in the latter 
to free rotators), the maximal Lyapunov exponent must vanish in these two 
limits. 

In the region of weak chaos, for U < 0.25, the curve has a weak A^-dependence. 
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Then Ai changes abruptly and a region of stronger chaos begins. In Ref. [1] 
it was observed that in between U = 0.2 and U = 0.3 a different dynamical 
regime sets in and particles start to evaporate from the main cluster. A sim- 
ilar regime was found in Ref. [9] and this behavior is also similar to the one 
found in Ref. [4] at a solid-liquid transition. In this region of strong chaoticity 
we observe a pronounced peak already for N — 100 [6], which persists and 
becomes broader for = 20000. 

The standard deviation of the kinetic energy fluctuations is plotted in 
Fig. 5(b) and it compares quite well with the theoretical prediction (15), al- 
though finite N effects are larger than for averaged quantities. In Fig. 5(c) we 
report the microcanonical specific heat obtained with formula (16), compared 
with numerical simulations at increasing values of N. Both these quantities 
display a similar behavior to the maximal Lyapunov exponent: they increase 
up to the critical point and then drop to zero. It is therefore quite natural to 
associate the growth of the Lyapunov exponent to the growth of fluctuations, 
as expressed both by the kinetic energy fluctuations and the specific heat. In 
this respect a similar connection was proposed in Ref. [17]. Unfortunately, a 
theoretical formula for the Lyapunov exponent as A" ^ oo does not yet exist 
(see anyway next section); however, the data in Fig. 5(a) already show what 
we should expect for the convergence of Ai as A^ increases. For [/ > ?7c, Ai — > 
as A^ — ^ cxD, thus revealing the presence of a whole region of integrability in this 
limit; rotators decouple, reducing the system to a "gas" of free rotators (which 
is consistent also with the vanishing of kinetic energy fluctuations). Fig. 6(a) 
shows that the convergence to zero can by fitted as Ai ~ N~^^^. This scaling 
law can be derived theoretically using a random matrix approximation [2] (see 
also next section), which is also shown in the same figure to approximate quite 
well the numerical results at large enough energy. (A N''^^^ convergence to the 
asymptotic value of the Lyapunov exponent is also observed in systems of hard 
spheres [34]). The question is still open whether Ai will show a discontinuity 
at Uc in the N ^ oo limit, as the kinetic energy fluctuations and speciflc heat 
do. Strong flnite size effects are present in the region U G [0.2,0.75]. This is 
the region where the cluster drifts (see [1]), while particles evaporate from and 
condensate on it. The Lyapunov exponent is here signiflcantly larger than at 
smaller energies. For U < 0.2, we observe a fast convergence to the A^ — > oo 
limit, and the scaling law Ai ~ VU is numerically obtained (see Fig. 6(b)). 
A heuristic justiflcation of this scaling was proposed in Ref. [2], and a new 
derivation is presented in Sec. 7.1. 

The extensivity of the Lyapunov spectrum was first proposed by Ruelle [35] 
and numerically tested in Refs. [36] (see [37] for a review). It amounts to 

check that plotting Aj vs. i/N and letting A^ go to infinity, while keeping fixed 
physically intensive parameters (in our case energy density or temperature), 
one obtains a convergence to an asymptotic curve, so called distribution of 
Lyapunov exponent. This is verified for our model in Fig. 7(a) at various energy 



11 



densities. The asymptotic Lyapunov distribution is more quickly reached, as 
N increases, at smaller energies (confirming the fast N convergence observed 
above for the maximal Lyapunov exponent at small energy). The data at 
U — 0.7 show a much weaker convergence to the asymptotic spectrum, as also 
confirmed by looking at the Kolmogorov- Sinai entropy density at the same 
value of U (upper points in the inset of Fig. 8). Heavier numerical simulations 
are needed to assess the convergence of the spectrum in this energy region. 
The spectra have a curious exponential shape in the first part (see the inset 
of Fig. 7(a)), which was already noticed in Ref. [18]. No significative change 
in the shape of the spectra at N — 50 and N — 100 is observed when going 
from below to above the phase transition point (see Fig. 7(b) for the N = 100 
case). However, we cannot exclude that differences in the spectra could arise 
at larger values of N, which are so far inaccessible to numerical experiments. 

The Kolmogorov- Sinai (K-S) entropy density (which, due to Pesin's formula, 
is in our case the sum of the positive Lyapunov exponents divided by N), 
Sks/N, is plotted in Fig. 8 against U. It shows again a peak at Uc, a fast 
convergence to a limiting value as N increases in the small energy limit and 
a slow convergence to zero (as expected) for U > Uc- The scaling laws are 
here roughly: Sks/N ~ C/^/^ at smaU U, see Fig. 9(a), and Sks/N ~ N'^/^ 
for U > Uc, see Fig. 9(b) (but a more refined numerical analysis is needed to 
confirm these results). 

Therefore, the analysis of the behavior of the Lyapunov spectrum and of the K- 
S entropy confirms the picture already emerging from the study of the maximal 

Lyapunov exponent, showing an increase of microscopic chaos near the phase 
transition point. It is somewhat intriguing that it is precisely when chaos is 
stronger that one reveals a slowing down of the relaxation to equilibrium (see 
Sec. 5); although not unexpected, because it is quite symplistic to associate 
the notion of chaos to that of efficient diffusion of orbits in phase space. When 
starting from a "water bag" initial condition one reaches the almost frozen 
QSS state, and also for this state the transient Lyapunov exponent has been 
checked to be positive, although slightly different than the asymptotic value 
in the equilibrium state. It could be that the phase space near the transition 
point has a rich structure with many coexisting chaotic QSS states. 



7 Analytical estimates of the mciximal Lyapunov exponent 

In this section we discuss some theoretical approaches which allow to justify 
the scaling laws of the maximal Lyapunov exponent observed in numerical 
experiments, and also to obtain some order of magnitude estimates. 
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7.1 Using the Vlasov solution 



One can try to use the stationary Vlasov distribution function (20,22) to derive 
some properties of the tangent space. For instance, averaging the Jacobian Jij 
(6,7) over this distribution, one obtains a constant matrix whose diagonal and 
off-diagonal elements are —M'^ + l/N and M'^/N, respectively (we are here also 
neglecting correlations among the particles). It is then quite easy, being now 
the Jacobian a constant matrix, to compute the maximal Lyapunov exponent. 
The result is 



with Lyapunov eigenvector {{S9i}, {Spi}) = (a, b), the vectors a and b being 
constant. This formula has the correct dependence on M, in fact at small 
energies the virial relation < K >r^< V > holds and then U = E/N ~ 
1 — M^, giving Ai ~ \/U , which is what is observed numerically. However, 
formula (28) predicts that the Lyapunov exponent vanishes with N^^^'^, which 
is not observed in numerical experiments. We have found numerically that, 
even if we allow for the true temporal fluctuations of M^, we obtain for the 
Lyapunov exponent the same value of formula (28) and the same Lyapunov 
eigenvector. On the contrary, if we look at the Lyapunov eigenvector given 
by the true dynamics, we observe that it is far from being constant, only 
a few components being significantly different from zero. We have checked 
that the number of nonvanishing components of the Lyapunov eigenvector 
remains constant as N increases (this has also been recently found also for 
a generalization of the HMF studied in Rcf. [23]). Thus, the typical size of 
each component of the vector remains constant as grows (remember that 
eigenvectors are normalized), while for constant eigenvectors this size decreases 
as N~^/'^] this is a naive argument that can explain the extra N~^/'^ present 
in formula (28) with respect to what is observed for the true Lyapunov, i.e. 
Ai ~ y/U. We have recently developed a perturbation theory scheme, which 
shows that for sparse eigenvectors Ai does not vanish with at small energies 
and, with some additional hypotheses, one obtains the y/U law [38]. 

7.2 A numerical test of a recently proposed formula 

We compare in this subsection a recently proposed formula for the maximal 
Lyapunov exponent [7,19] with numerical results, giving also an intuitive in- 
terpretation of it for our model. In Ref. [7] a general method is formulated to 
describe Hamiltonian chaos using the differential geometric structure underly- 
ing the dynamics. They consider the Hamiltonian many-body dynamics as a 




(28) 
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geodesic flow on a Riemannian manifold; then the chaotic motion reflects the 
instabihties of the geodesic flow, which depend on the curvature properties of 
the manifold. They obtain the following formula for Ai 

1 

= • A=(2<'^- + \/(f^)'+P-aT)j' , (29) 



where Qq ^-nd are the average Ricci curvature and the variance of its fluctu- 
ations respectively. For the HMF, these two quantities turn out to be given by 
very simple expressions in terms of the averages and fluctuations introduced 
in Sec. 3 (see also ref. [19]) 



1]o = <M2>-1 = T, + (1-2C/)-1 

a'n^Nal,,^4El, (30) 

where (7^2 is the variance of the fluctuations of M^. All this is the consequence 
of the fact that the Ricci curvature,— 1/iV J2i Ju with Ja given by formula (6), 
is simply given in our model by — 1/N. The Ricci curvature is assumed 
in Ref. [7] to be a 5-correlated stochastic process in time, and the correlation 
time T in formula (29) is estimated as 

r- = 2(rr^ + r.-); r, ^ \ ; r. = ^ (31) 



Some comments about the derivation and the use of formula (29) are very 
important. It is derived in the "diagonal" approximation, which corresponds 
to neglect the effect of off-diagonal terms (7) in Eqs. (5). This approximation 
has been checked numerically not to be valid at small energies, leading to 
a value of Ai which is orders of magnitude less than the true one (see also 
Fig. 10). The assumption of a 5-correlated stochastic process makes easier the 
calculation of Ai, but reduces the range of applicability of formula (29). The 
estimate for r is a rather delicate problem, where some arbitrariness can enter 
the theory [39]. 

In the high energy phase, M fluctuates above zero and scales with A^~^/^ at 
large A^, then we have < >~ (7^2 ~ -/V~^ and moreover r ~ T2/2 since 
T2 = M/{aM^VN) = 0(1) while n ~ A^^/l Hence one flnds [19] 

Ai ~ A ~ {4:Nal^2Tf'^ ~ Ar-V3 _ ^33) 



14 



In the low energy phase (7m2 <^ and Ai reduces to 



^1 ~ ^ (33) 



which gives a relation betwen Ai and the fluctuations of kinetic energy. This 
result supports the link between chaos and kinetic energy fluctuations already 
claimed in Ref.[2] although the formula derived heuristically there was differ- 
ent, but in better agreement with the scaling law for small U. 

In Fig. 10 we compare the numerically computed values of Ai at increasing 
values of with those obtained from formula (29), using both the finite N 
numerical values for < > and and the thermodynamic limit values 
in the microcanonical ensemble. We have also checked that cqs. (30) are well 
reproduced by numerical simulations. Formula (29) reproduces the behavior 
of Ai vs. U within a factor of two over a large range of energies. In particular 
it predicts correctly a maximum of Ai around Uc and the A'"~^/^-law in the 
high energy phase. 

However, in the limit ^ this formula gives Ai ~ C/^, as can be easily 
derived from formula (33), which is sharply in contrast with the behaviour 
Ai ~ \/U observed in numerical simulations (sec Fig. 6(b)). We think that, in 
this limit, the stochastic approximation for the average Ricci curvature breaks 
down, because this quantity is in our case nothing but M^, and this is an al- 
most regularly oscillating quantity as C/ — > 0, due to the collective excitations 
of the cluster. However, it is possible to use r as a fitting parameter to repro- 
duce the numerical data. In fact some preliminary numerical investigations 
have shown that in this region the correlation time given by eq. (31) strongly 
underestimates the realistic one. 



8 CONCLUSIONS 



We have investigated the dynamical and statistical behavior of a system with 
long-range forces showing a second order phase transition. Both the maximal 
Lyapunov exponent Ai and the Kolmogorov-Sinai entropy density Sks/N are 
peaked at the phase transition point, where kinetic energy fluctuations and 
specific heat are maximal. There is actually a small shift to lower energies 
due to finite size effects. The latter are present also in the Lyapunov spec- 
tra and in the Kolmogorov-Sinai entropy. Above the phase transition point, 
both Ai and Sks vanish as iV ^ oo. We think that this toy model contains 
some important ingredients to understand the behavior of macroscopic order 
parameters when dynamical chaos is present at the microscopic level. Most 
of our findings are probably common to other Hamiltonian systems showing 
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second order phase transitions. In particular our results could be very impor- 
tant in order to understand the relaxation to the equilibrium solution and the 
success of statistical approaches in describing the nuclear multifragmentation 
phase transition. 
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Fig. 1. Theoretical predictions in the canonical ensemble (full curve) for the mag- 
netization M vs. U, panel (a), and for the caloric curve T vs. U, panel (b), in 
comparison with numerical simulations (microcanonical ensemble) for N= 100, 1000, 
5000, 20000. The vertical line indicates the critical energy Uc = 3/4. We plot also the 
microcanonical results for the Quasi-Stationary States (QSS) in the case N=20000 
(losanges) in the inset of panel (b). 

Fig. 2. Equilibrium reduced distribution functions in angles 9 and momenta p for 
U = 0.095 ((a) and (b)) and U = 0.36 ((c) and (d)). The numerical simulations for 
= 1000 (histogram) are compared with the theoretical distributions (20,22) 

Fig. 3. Reduced distribution function in p fov U = 0.69 and N = 1000 at increasing 
times. The initial distribution is a water bag (histogram). The thick full line is the 
equilibrium distribution. 

Fig. 4. (a) Relaxation to equilibrium of the Boltzmann entropy of the reduced 
distribution in p for N=1000 and N=5000; (b) Linear increase with N of the time 
scale Tr for the relaxation to equilibrium 

Fig. 5. (a)Numerical data for the largest Lyapunov exponent as a function of U for 
various system sizes: N=100,1000, 5000 and 20000. (b) Kinetic energy fluctuations 
vs. U. (c) Specific heat obtained from formula (16). The vertical line indicates the 
critical energy Uc = 3/4. The full lines in panels (b) and (c) are the microcanonical 
predictions. 

Fig. 6. (a) Convergence to zero of the largest Lyapunov exponent vs. N in the high 
energy region. The fitted curve is the power law N~^/^, the circles are the result of 
a random matrix simulation, (b) U^^'^ growth of the maximal Lyapunov exponent 
vs. U at small U, a rather weak N dependence is observed. 

Fig. 7. (a)Lyapunov spectra for three different values of U at different N in linear 
and lin-log scale (inset for N = 20) showing the convergence (weaker for the U = 0.7 

case) to the asymptotic distribution. (b)Lyapunov spectra for fixed N = 100 with 
two different values of the energy, below and above the critical energy. No significa- 
tive differences are observed at these values of N. 

Fig. 8. The Kolmogorov-Sinai entropy density Sks/^ vs. U at increasing values 
of N . The convergence to the thermodynamic limit of the K-S entropy density is 
shown in the inset for U = 0.01,0.1,0.3,0.7. The convergence is weaker for this latter 
U value. 

Fig. 9. (a) Growth of the K-S entropy density at small U with the exponent 3/4. A 

rather weak A^-dependence is observed, (b) Convergence to zero of the K-S entropy 

density as N increases in the high energy region. The fitted curve is the power law 
iV-V5. 
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Fig. 10. Comparison of formula (29) for Ai at the thermodynamic hmit (full line) and 
for finite N (shaded circles) with numerical simulations (open circles) at increasing 
values of N. 
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